/* ****************************************************************************
 *
 * Copyright (c) Microsoft Corporation. 
 *
 * This source code is subject to terms and conditions of the Microsoft Public License. A 
 * copy of the license can be found in the License.html file at the root of this distribution. If 
 * you cannot locate the  Microsoft Public License, please send an email to 
 * dlr@microsoft.com. By using this source code in any fashion, you are agreeing to be bound 
 * by the terms of the Microsoft Public License.
 *
 * You must not remove this notice, or any other, from this software.
 *
 * THIS CODE AND INFORMATION ARE PROVIDED "AS IS" WITHOUT WARRANTY OF ANY
 * KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A
 * PARTICULAR PURPOSE.
 * ***************************************************************************/


/*! \brief A Fast Fourrier Transform template (to use with float or double)

Do not use this class with other types than float or double, it will lead to unpredictable results

This class computes the FFT of a vector.

\par Credits: Original code

    <p>Don Cross <dcross@intersrv.com></p>

    <p>http://www.intersrv.com/~dcross/fft.html</p>

    <p>Contains definitions for doing Fourier transforms
    and inverse Fourier transforms.
    This module performs operations on arrays of 'float' and 'double'</p>

\par Modification
	<p>Encapsulated all in a C++ template class and added STL support...</p>

\ingroup AlgoToolsGroup
*/
template<class T>
class TFFT
{
public:

	/*! FFT Transform

	\param bInverseTransform  true=forward FFT, false=inverse FFT
	\param vRealIn array of input's real samples
	\param vImagIn array of input's imag samples 
	\param vRealOut array of output's reals 
	\param vImagOut  array of output's imaginaries 

	If vImagIn is empty, on vRealIn is transformed

	*/
	static bool Transform(bool bInverseTransform, 
	const std::vector<T>&   vRealIn,        
    const std::vector<T>&   vImagIn,   
    std::vector<T>&   vRealOut,   
    std::vector<T>&   vImagOut );   

	/*! FFT Transform,
	\param bInverseTransform  true=forward FFT, false=inverse FFT
	\param vRealIn array of input's real samples
	\param vImagIn array of input's imag samples 
	\param vRealOut array of output's reals 
	\param vImagOut  array of output's imaginaries 

	If vImagIn is empty, on vRealIn is transformed
	*/
	static bool Transform(bool bInverseTransform, 
	const std::valarray<T>&   vRealIn,        
    const std::valarray<T>&   vImagIn,   
    std::valarray<T>&   vRealOut,   
    std::valarray<T>&   vImagOut );   

	//! return true if x is a power of 2
	inline static bool IsPowerOfTwo( unsigned x );
	/*!   \ brief Abstract frequency
		The following function returns an "abstract frequency" of a
	given index into a buffer with a given number of frequency samples.
	Multiply return value by sampling rate to get frequency expressed in Hz.
	*/
	inline static double IndexToFrequency( unsigned NumSamples, unsigned Index );

private:
	inline static unsigned NumberOfBitsNeeded( unsigned PowerOfTwo );
	inline static unsigned ReverseBits( unsigned index, unsigned NumBits );

	static double DDC_PI; 
};

template<class T>
double TFFT<T>::DDC_PI=3.141592653;

template<class T>
inline bool TFFT<T>::IsPowerOfTwo(unsigned x )
{
    if ( x < 2 )
        return false;
    if ( x & (x-1) )        // Thanks to 'byang' for this cute trick!
        return false;
    return true;
}

template<class T>
inline unsigned TFFT<T>::NumberOfBitsNeeded ( unsigned PowerOfTwo )
{
    unsigned i;
	ASSERT( PowerOfTwo > 2 );
    for ( i=0; ; i++ )
    {
        if ( PowerOfTwo & (1 << i) )
            return i;
    }
}

template<class T>
inline unsigned TFFT<T>::ReverseBits( unsigned index, unsigned NumBits )
{
    unsigned i, rev;
    for ( i=rev=0; i < NumBits; i++ )
    {
        rev = (rev << 1) | (index & 1);
        index >>= 1;
    }

    return rev;
}

template<class T>
inline double TFFT<T>::IndexToFrequency( unsigned NumSamples, unsigned Index )
{
    if ( Index >= NumSamples )
        return 0.0;
    else if ( Index <= NumSamples/2 )
        return (double)Index / (double)NumSamples;

    return -(double)(NumSamples-Index) / (double)NumSamples;
}

template<class T>
bool TFFT<T>::Transform(bool bInverseTransform, 
	const std::vector<T>&   vRealIn,        
    const std::vector<T>&   vImagIn,   
    std::vector<T>&   vRealOut,   
    std::vector<T>&   vImagOut )  
{
    unsigned NumBits;    /* Number of bits needed to store indices */
    unsigned i, j, k, n;
    unsigned BlockSize, BlockEnd;
	unsigned NumSamples;
	
	if (vImagIn.empty())
		NumSamples=vRealIn.size();
	else
		NumSamples= __min( vRealIn.size(), vImagIn.size());

    double angle_numerator = 2.0 * DDC_PI;
    double tr, ti;     /* temp real, temp imaginary */

    if ( !IsPowerOfTwo(NumSamples) )
		return false;

	vRealOut.resize(NumSamples);
	vImagOut.resize(NumSamples);

    if ( bInverseTransform )
        angle_numerator = -angle_numerator;

    NumBits = NumberOfBitsNeeded ( NumSamples );

    /*
    **   Do simultaneous data copy and bit-reversal ordering into outputs...
    */

    for ( i=0; i < NumSamples; i++ )
    {
        j = ReverseBits ( i, NumBits );
        vRealOut[j] = vRealIn[i];
		vImagOut[j] = ( vImagIn.empty() ) ? (T)0.0 : vImagIn[i];
    }

    /*
    **   Do the FFT itself...
    */

    BlockEnd = 1;
    for ( BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1 )
    {
        double delta_angle = angle_numerator / (double)BlockSize;
        double sm2 = sin ( -2 * delta_angle );
        double sm1 = sin ( -delta_angle );
        double cm2 = cos ( -2 * delta_angle );
        double cm1 = cos ( -delta_angle );
        double w = 2 * cm1;
        double ar[3], ai[3];
        //double temp;

        for ( i=0; i < NumSamples; i += BlockSize )
        {
            ar[2] = cm2;
            ar[1] = cm1;

            ai[2] = sm2;
            ai[1] = sm1;

            for ( j=i, n=0; n < BlockEnd; j++, n++ )
            {
                ar[0] = w*ar[1] - ar[2];
                ar[2] = ar[1];
                ar[1] = ar[0];

                ai[0] = w*ai[1] - ai[2];
                ai[2] = ai[1];
                ai[1] = ai[0];

                k = j + BlockEnd;
                tr = ar[0]*vRealOut[k] - ai[0]*vImagOut[k];
                ti = ar[0]*vImagOut[k] + ai[0]*vRealOut[k];

                vRealOut[k] = vRealOut[j] - (T)tr;
                vImagOut[k] = vImagOut[j] - (T)ti;

                vRealOut[j] += (T)tr;
                vImagOut[j] += (T)ti;
            }
        }

        BlockEnd = BlockSize;
    }

    /*
    **   Need to normalize if inverse transform...
    */

    if ( bInverseTransform )
    {
        double denom = (double)NumSamples;

        for ( i=0; i < NumSamples; i++ )
        {
            vRealOut[i] /= (T)denom;
            vImagOut[i] /= (T)denom;
        }
    }

	return true;
}
